function hypcompare(varargin)
%
% compare hyp funcitons from the dynamic clamp experiment.
% Plots the FSLs from a group of experiments
% Paul B. Manis, Ph.D.
% July 25-27, 2007. 
% for Street's paper on rHpTx, dynamic clamp, immuno and insitu.
%

if(nargin == 0)
    cellid = 1;
else
    cellid = varargin{1};
    if(isempty(cellid))
    cellid = [1,2,3,5,6,7]; % get the ones used in the paper
    end;
end;
if(nargin > 1)
    targetv = varargin{2};
else
    targetv = -70;
end;
if(nargin > 2)
    makefileflag = 1;
    fprintf(1, 'Making data files from RTXI text files\n');
else
    makefileflag = 0;
end;

if(makefileflag == 1 )
    mon = {'jan', 'feb', 'mar', 'apr', 'may', 'jun', 'jul', 'aug', 'sep', 'oct', 'nov', 'dec'};
    for j = 1:length(cellid)
        [cellname, leg, explist] = setcellinfo(j);
        for i = 1:length(explist)
            fn = sprintf('%s/exp%03d/Hyp_Curve', cellname, explist{i});
            % convert cellname to our standard
            [yr res] = strtok(cellname, '.');
            [mo res] = strtok(res, '.');
            [day res] = strtok(res, '.');
            if(isempty(res))
                cn = 'z';
            else
               [ cn, res] = strtok(res, '.');
            end;

[s bp0] = system('echo $HOME');
basep = strcat(bp0, '/Desktop/Street-rHpTx-DC/');
            fnout = sprintf('%02d%3s%02d%c_%03d.pbm', str2num(day), mon{str2num(mo)}, mod(str2num(yr), 100), cn, explist{i});
%             fnout
%             [basep fnout]
            rtxi2datac([basep fn], [basep fnout], [leg{i}]);
        end;
    end;
    return;

end;

plotflag = 0;
summary=[];

for i = 1:length(cellid)
    [hfig, summary] = hypcompare2(cellid(i), targetv, summary);
    if(plotflag)
        orient landscape;
        print (hfig, '-dpsc2');
    end;
end;
% now generate a summary of the relative shifts.
% express
if(length(cellid) > 2)
    hn = newfigure('summary', sprintf('summary of %d cells', length(cellid)));

    for i = 1:length(cellid)
        rmp(i,:) = summary(i).rmp;
        fsl(i,:) = summary(i).fsllatv;
        rmpd(i,:) = summary(i).rmp - summary(i).rmp(1);
        fsld(i,:) = summary(i).fsllatv - summary(i).fsllatv(1);
        fslp(i,:) = 100*summary(i).fsllatv/summary(i).fsllatv(1);
        mrd(i,:) = summary(i).mr - summary(i).mr(1);
        mrp(i,:) = 100*summary(i).mr/summary(i).mr(1);
        mr(i,:) = summary(i).mr;
    end;

    mrmpd = mean(rmpd, 1);
    mfsld = mean(fsld, 1);
    mfslp = mean(fslp, 1);
    mmrd = mean(mrd, 1);
    mmrp = mean(mrp, 1);
    srmpd = std(rmpd, 1);
    sfsld = std(fsld, 1);
    sfslp = std(fslp, 1);
    smrd = std(mrd, 1);
    smrp = std(mrp, 1);
    srmpd(1) = std(rmp(:,1));
    sfsld(1) = std(fsl(:,1));
    sfslp(1) = 100*sfsld(1)/mfsld(1);
    smrd(1) = std(mr, 1);
    smrp(1) = 100*smrd(1)/mmrd(1); % convert to percentage

    xl = [1:size(rmpd, 2)];
    subplot(2,3,1);
    errorbar(xl, mrmpd, -srmpd, srmpd, 'kx-');
    subplot(2,3,2);
    errorbar(xl, mfsld, -sfsld, sfsld, 'ko-');
    subplot(2,3,5);
    errorbar(xl, mfslp, -sfslp, sfslp, 'ro-');
    subplot(2,3,3);
    errorbar(xl, mmrd, -smrd, smrd, 'ks-');
    subplot(2,3,6);
    errorbar(xl, mmrp, -smrp, smrp, 'rs-');
    disp 'condition'
    xl'
    disp 'mean rmp'
    [mrmpd; srmpd]'
    disp 'fsl delta'
    [mfsld; sfsld]'
    disp 'fsl percentile'
    [mfslp; sfslp]'
    disp 'isi delta'
    [mmrd; smrd]'
    disp 'isi percentile'
    [mmrp; smrp]'

    % dump the raw result matrices
    disp 'fslatv'
    [summary.fsllatv]'
    disp 'meanrate'
    [summary.mr]'
    disp 'rmp'
    [summary.rmp]'

end;


function [hf, summary] = hypcompare2(cellid, targetv, summary, bpath)

global t v


[cellname, leg, explist] = setcellinfo(cellid);

symlist='os^dvph><x.*+';
nsym = length(symlist);
this_sym = mod([0:length(explist)-1], nsym)+1;
basep = [bpath cellname '/'];
hf = newfigure('hypcompare', 'Comparision of Hyp protocols');
figure(hf);
clf;
hm = multiplot(hf, 2, 3, 1, [0 0 0.5 1], [2 5]); % put all the plots on the left side of the figure

%leg = {};
misi=zeros(length(explist), 1);
sisi = misi;
cond = misi;
tstep = 110;
mr = zeros(length(explist), 1);
sr = zeros(length(explist), 1); % preallocate
trlist = zeros(length(explist), 1);
rmp = zeros(length(explist), 1);
title(basep);
vlen = 0;
for i = 1:length(explist)
    v=[]; t=[];

    fn = sprintf('exp%03d/Hyp_Curve', explist{i});
    leg2{i} = sprintf('exp%03d  ', explist{i});
    [v,curr,t,d] = rtxi_read([basep, fn]);
    if(isempty(d))
        return;
    end;
    v = v * 1000; % scale to mV
    t = t * 1000; % scale to msec
    if (i == 1)
        vlen = size(v,2);
        vlen = 9200;
        vtr = zeros(length(explist), vlen); % preallocate room for traces
        ttr = vtr;
    end;

    vpre = getvalue('mean', v, d, 80, 100); % convert to mV
    rmp(i) = mean(getvalue('mean', v, d, 0, 2));

    [itr, jtr] = min(abs(vpre-targetv));
    trlist(i) = jtr;
    vtr(i, :) = v(jtr, 1:vlen);
    ttr(i, :) = t(jtr, 1:vlen);

    %     figure
    %     plot(t', v');
    %     drawnow

    % get fsl data
    [fsl, fisi, nsp, spikes, plkat] = find_spikes2(d, v, 0, 1000, -20);

    % organize fsl data
    dx=[vpre;fsl-tstep]';
    dy=sortrows(dx, 1);
    figure(hf);
    axes(hm(1)); % 1 is the FSL plot

    mysym = sprintf('%c-', symlist(this_sym(i)));
    plot(dy(:,1), dy(:,2), mysym, 'markersize', 8);
    set(gca, 'Ylim', [0 200]);
    hold all
    xlabel('Vpre (mV)');
    ylabel('FSL (ms)');

    fslatv(i) = dy(jtr, 2); % save fsl measured at the selected V

    % compute isis
    misi=zeros(length(spikes),1);
    sisi=zeros(length(spikes),1);
    cond=zeros(length(spikes),1);
    for j = 1:length(spikes)
        misi(j) = mean(diff(spikes(j).latency));
        sisi(j) = std(diff(spikes(j).latency));
        cond(i) = i;
    end;
    %figure(hf2);
    axes(hm(2)); % 2,3,2 is the mean ISI plot
    he=errorbar(dy([spikes.source],1), misi', sisi, -sisi, mysym);
    hold all
    xlabel('Vpre (mV)');
    ylabel('Mean ISI (ms)');
    set(gca, 'Ylim', [0 100]);
    mr(i) = mean(misi');
    sr(i) = std(misi); % standard deviation across all voltages
end;
figure (hf);
axes(hm(1)); % make a legend for the fsl plot, but put it away from plot
d=cell(size(leg,2), 1);
d{1} = sprintf('Cell: %s', cellname);
for i = 1:size(leg, 2)
    d{i} = sprintf('%s:%s', char(leg2(i)), char(leg(i)));
end;
legend(d);


% info text box.
axes(hm(6)); % subplot('Position', [0.70, 0.8, 0.24, 0.15]);
axis([0,1,0,1])
axis('off')
text(0, 0.8, sprintf('Cell: %s ', cellname));
% now move the legend to a clear area
hl =findobj('tag', 'legend');
h6pos = get(hm(6), 'Position');
h6pos(2) = h6pos(2)*1.1;
h6pos(4) = h6pos(4)*2/3;
set(hl, 'Position', h6pos);


bx = [1:length(mr)];

axes(hm(3)); % fsl versus condition here
bar(bx, fslatv);
hold on;
%plot(bx, fslatv, '-s');
text(bx, fslatv, leg, 'horizontalalignment', 'left', 'rotation', 65,...
    'verticalalignment', 'middle', 'margin', 10, 'fontsize', 10);
ylabel(sprintf('FSL (ms) at %2d mV', targetv));
xlabel('Condition');

axes(hm(4)); % mean isi plot goes here
bar(bx, mr);
hold on;
errorbar(bx, mr, sr, -sr, '-');
text(bx, mr, leg, 'horizontalalignment', 'left', 'rotation', 65,...
    'verticalalignment', 'middle', 'margin', 10, 'fontsize', 10);
% set(gca, 'XtickLabel', leg);
ylabel('Mean ISI (ms)');
xlabel('Condition');


% make legend box
d=cell(size(leg,2), 1);
for i = 1:size(leg, 2)
    d{i} = sprintf('%s:   %s', char(leg2(i)), char(leg(i)));
end;
%ha=annotation('textbox', [0.5, 0.1, 0.45, 0.4]);
%set(ha, 'String', d);

% subplot(2,3,6);  % RMP plot
% bar(bx, rmp);
% hold on;
% %plot(bx, rmp, '-s');
% % text(bx, rmp, leg, 'horizontalalignment', 'left', 'rotation', 65,...
% %     'verticalalignment', 'middle', 'margin', 10, 'fontsize', 10);
% ylabel('RMP (mV)');
% xlabel('Condition');

axes(hm(5));  % RMP plot
hold on;
co = get(gca, 'colororder');
for i = 1:length(explist)
    mysym = sprintf('-%c', symlist(this_sym(i)));
    plot(rmp(i), fslatv(i), mysym, 'markersize', 8, 'markeredgecolor', co(i,:), ...
        'markerfacecolor', 'none')
    plot(rmp(i),  mr(i), mysym, 'markersize', 8, 'markeredgecolor', co(i,:), ...
        'markerfacecolor', co(i,:))
end;
plot(rmp, fslatv, 'k-');
plot(rmp, mr, 'b-');
% text(bx, rmp, leg, 'horizontalalignment', 'left', 'rotation', 65,...
%     'verticalalignment', 'middle', 'margin', 10, 'fontsize', 10);
xlabel('RMP (mV)');
ylabel('fsl or fisi (ms)');
% now plot the selected traces for the target voltage
stackplot(hf, ttr, vtr, 1, [0.52 0 0.47 1]);
%plot(ttr', vtr');

%
% generate some summary information for this cell
n = length(summary)+1;
summary(n).rmp = rmp;
summary(n).fsllatv = fslatv;
summary(n).mr = mr;
summary(n).targetv = targetv;


%----------------------
% set the information about this particular selection;;;

function [cellname, leg, explist] = setcellinfo(cellid)

switch(cellid)
    case 1
        cellname = '2007.4.25'; % ok
        leg={'control', 'toxin', 'IKIF(800)', 'IKIF(800)+IKIS(40)', 'IKIS(30)', 'no DC'};
        explist = {2 11 12 14 15 16};

    case 2
        cellname = '2007.4.27.a'; % possible exemplar cell...
        % 2007.4.27a compare: control, toxin, IKIF(1000), IKIF(100)+IKIS(30), IKIS(30) no DC
        leg=({'control', 'toxin', 'IKIF(1000)', 'IKIF(1000)+IKIS(43)', 'IKIS(30)', 'no DC'});
        explist = {4 16 20 17 22 23};
    case 3
        cellname = '2007.4.27.b'; % toxin action not clear at -75, ok at -65. Used as example in paper
        leg=({'control', 'toxin', 'IKIF(800)', 'IKIF(800)+IKIS(30)', 'IKIS(30)', 'no DC'});
        explist = {4 6 7 11 10 13};
    case 4 % not good cell - not much fsl shift after toxin on. Little effect of dc
        cellname = '2007.5.8.a';
        leg=({'control', 'toxin', 'IKIF(1000)', 'IKIF(1000)+IKIS(175)', 'IKIS(175)', 'no DC'});
        explist = {12 15 17 19 20 24};
    case 5
        cellname = '2007.5.17.a'; % ok, small shifts, clear pattern
        leg=({'control', 'toxin', 'IKIF(700)', 'IKIF(700)+IKIS(30)', 'IKIS(30)', 'no DC'});
        explist = {2 6 8 10 11 12};

    case 6 % second cell in file group from 2007.5.17.a. Files moved in directory to b (see pdf file of scanned notes)
        cellname = '2007.5.17.b'; % ok, small shifts, clear pattern
        leg=({'control', 'toxin', 'IKIF(700)', 'IKIF(700)+IKIS(30)', 'IKIS(30)', 'no DC'});
        explist = {20 24 25 26 27 29};

    case 7
        cellname = '2007.5.25'; % ok, small shifts
        leg=({'control', 'toxin', 'IKIF(500)', 'IKIF(500)+IKIS(30)', 'IKIS(30)', 'no DC'});
        explist = {2 5 6 7 8 9};
    case 8
        cellname = '2007.5.25.b'; % ok, good shifts, but cell has ADP in AHP- CW?.
        leg=({'control', 'toxin', 'IKIF(400)', 'IKIF(400)+IKIS(30)', 'IKIS(30)', 'no DC'});
        explist = {2 5 6 7 8 9};
    otherwise
        return;
end;
